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Abstract. In this paper we improve traditional steepest descent methods for the direct min- 
imization of the Gross-Pitaevskii (GP) energy with rotation at two levels. We first define a new 
inner product to equip the Sobolev space and derive the corresponding gradient. Secondly, for 
the treatment of the mass conservation constraint, we use a projection method that avoids more 
complicated approaches based on modified energy functionals or traditional normalization methods. 
The descent method with these two new ingredients is studied theoretically in a Hilbort space setting 
and wc give a proof of the global existence and convergence in the asymptotic limit to a minimizer 
of the GP energy. The new method is implemented in both finite difference and finite element two- 
dimensional settings and used to compute various complex configurations with vortices of rotating 
Bosc-Einstein condensates. The new Sobolev gradient method shows better numerical performances 
compared to classical or gradient methods, especially when high rotation rates arc considered. 
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1. Introduction. First experimental realizations of Bose-Einstein condensates 
(BECs) in dilute alkali-metal gases [U [121 HI] led to an explosion of mathematical 
and theoretical studies aimed at better understanding such systems. Recent efforts 
were devoted to documenting the superfluid nature of the condensate by providing 
evidence for the existence of quantized vortices when rotating the condensate. It was 
indeed experimentally observed [25l [26l [U [19] that instead of solid body rotation, the 
condensate rotates by forming vortices with quantized circulation. Initially, a few 
vortices are formed, and, as the rotation frequency increases, the vortices form an 
array similar to the Abrikosov lattice observed in type II superconductors. Since the 
rotating BEC is a highly controllable system with a simple theoretical description, 
it provides a perfect set-up for the theoretical study of macroscopic systems with 
quantized vortices. 

In the zero-temperature limit, a dilute gaseous BEC is mathematically described 
by a macroscopic wave function derived in the framework of the Gross-Pitaevskii 
(GP) mean field theory. The spatial configuration of the wave function ipi'^)^ with 
X = {x, y, zY , is obtained by minimizing the GP energy in the rotating frame, 

;^|VV|' + V^trap|^|' + ||^|' + iftrrj-(xx V)^, (1.1) 
Jr3 2m 2 

subject to the normalization condition, \ip\'^ = N, with N the number of particles 
(atoms). In the previous expression, h is Planck's constant, m the atomic mass of 
the gas, n the angular velocity vector, and Vtrap the magnetic trapping potential 
with trap frequencies {ujx,ujy,ujz). Wc denote by ip* the complex conjugate of ip. 
The interactions between atoms are described hy g ~ , with Ug the s-wave 

scattering length. As in most of experimental settings, we consider in the following 
that fl = (0, 0, and that Vtrap has a lower bound and Vtrap(x) oo, as x — >• cxd. 
Since from the previous assumption we can infer that 'ip{'X-) ^ 0, as x — > oo, it 
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suffices to work in a bounded domain T) C with homogeneous Dirichlet boundary 
conditions = on dD. 

In practice it is common to scale the energy so that the units become dimension- 
less. Using the scaUng r = x/d, w(r) = ^{k)cP/^ /^/N, fl ~ fl/ujj_, with d = 



the harmonic-oscillator length and uj± = min(a;a;, i^y) the transverse trap frequency, 
we obtain the non-dimensional energy (per particle) functional: 

Eiu) = / + Vtrap\u\^ + - niu*iA'V)u, (1.2) 

where Vtrap = n^^trap, 9 = , and A = {y,—x,0). The mass conservation 

constraint becomes: 



u\' = \\ur = i, (1.3) 

where we denote by ||.|| = || .llis^p^^). 

For given constants 51,(7, S'Hd trapping potential function Vtrap, the minimizer 
Ug of the functional (jl.2p under the constraint (|1.3p is called the ground state of the 
condensate. Local minima of the energy functional with energies larger that Eiug) are 
called excited (or metastable) states of the condensate. For a detailed discussion of 
the derivation of the Gross-Pitaevskii energy and the physics of rotating Bose-Einstein 
condensates see, for example, [15] and [24] . 

The two key issues in numerically computing ground or excited states of BEG 
are (i) how to derive a numerical algorithm that starts from a chosen initial state 
and iteratively diminishes the energy of the solution to rapidly converge to a local 
minimum of the functional (jl.2[) . and (ii) how to take into account the mass constraint 
(|1.3p . These two issues are obviously connected and have to be considered together 
in deriving efficient numerical algorithms. We present in this paper new approaches 
to address both issues and prove their superior numerical performance in the case of 
the energy minimization of the Gross-Pitaevskii energy with rotation. 

Most of the numerical algorithms proposed in the literature use the so-called 
normalized gradient flow [9], that consists in two steps: the steepest descent method 
is applied to the unconstrained problem, 

du ldE{u) V^u , ,2 , .^,4, 



dt 2 du 



VtrapU- g\u\ u + i9.A Vu, (1.4) 



to advance the solution from the discrete time level tn to tn+i; the obtained predictor 
'!i(r, i„+i) is then normalized in order to satisfy the unitary norm constraint and set 
the solution at i„+i: 

M(r,i„+i = — — (1.5 

l|M(r,t„+i)|| 

The gradient flow equation (|1.4p (or the related continuous gradient flow equation, 
see [9]) can be viewed as a complex heat equation and, consequently, solved by different 
classical time integration schemes (Runge-Kutta-Fehlberg [17] , backward Euler [9l [S] 
[To] , second-order Strang time-splitting [9l[5], combined Runge-Kutta-Grank-Nicolson 
scheme [3l|4l[T3], etc.), and different spatial discretization methods (Fourier spectral 
[T7] . finite elements [5], finite differences [HI [31 HI [13], sine-spectral [5], Laguerre- 
Hermite pseudo-spectral [TU], etc.). 
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It is interesting to note that in the descent method (jl.4p . the right-hand side 
represents the L^-gradient (or ordinary gradient) of the energy functional An impor- 
tant improvement of the convergence rate of the descent method was obtained in [T7] 
by replacing the ordinary gradient with the gradient defined on the Sobolev space 
i?^(r>, C). The same Sobolev gradient method (see [27] for various applications of 
this method) was recently used to minimize simpler Schrodinger type functionals in 
[29j . Similar increase of the convergence rate over the ordinary gradient method was 
reported. A first new contribution of the present paper is to introduce a new definition 
of the inner product to equip the Sobolev space in the case of the GP energy with 
rotation. A proof of the existence of the asymptotic limit for the evolution equations 
associated with the Sobolev gradients in a Hilbert space setting is also given. When 
implemented in a finite difference or finite element settings, the new Sobolev gradient 
method shows better numerical performances compared to classical or gradient 
methods, especially when high rotation rates (f2) are considered. 

The second important contribution of this work concerns the issue of the mass 
conservation constraint (jl.3p . Instead of the classical (and very popular) normaliza- 
tion approach (|1.5p . we suggest a projection method that preserves the norm of the 
initial state through the minimization procedure. The idea to project the (Sobolev) 
gradient into the tangent space associated to the constraint was already used to de- 
rive numerical algorithms for minimizing harmonic maps [SJ I28| . and, recently, to 
numerically find the smallest eigenvalue and corresponding eigenvectors of a Hermi- 
tian operator [7] . Different algorithms based on the projected gradient were developed 
in these studies and successfully applied to different energy functionals: the Oseen- 
Frank energy for liquid crystals |6] , the Dirichlet energy of harmonic maps |28| and the 
Hartree-Fock energy for quantum chemical molecular systems [7|. We derive here a 
projected Sobolev gradient method adapted to the Gross-Pitaevskii energy functional 
and provide an explicit expression of the projected gradient that allows to minimize 
trajectories when Hilbert spaces other than are considered. The new projection 
method proved very helpful in numerical implementations and allowed to avoid al- 
ternative methods to treat the mass constraint by adding to the energy functional a 
penalty term with a Lagrange multiplier (e. g. |17[ IllQ . 

The organization of the paper is as follows. In section[2]we introduce an alternate 
inner product on the Sobolev space and show that this inner product is equivalent 
to the traditional inner product on H^. The corresponding new Sobolev gradient 
is also derived. We discuss in section |3] a constructive projection method for the 
mass constraint and give our existence and convergence result for the asymptotic 
limit of the evolution equation defined by the Sobolev gradients. In section [¥] we 
give a discussion of the finite difference and finite element implementations in two- 
dimensions. The last section is devoted to numerical tests designed as benchmarks 
to compare performances of different Sobolev gradient methods. The effectiveness 
of the newly proposed Sobolev gradient method is proved by computing stationary 
states of rotating BEG that are physically relevant (high rotation and large interaction 
constants). 

2. Gradient descent methods using several gradients. In optimization 
problems that use a gradient descent or ascent technique, one usually has a choice of 
norms to use in the argument. If the norm has an associated inner product, then one 
can obtain a gradient with respect to this inner product (see |27j for an explanation). 
In the example of the minimization problem of Schrodinger type functionals, the 
gradient represents the direction of change per unit time. Therefore, one wants to 
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choose the gradient in the descent method (|1.4[) so that the change in energy is 
maximal at each step. For the case of the Gross-Pitaevskii energy with rotation, we 
notice that the energy (|1.2p can be written as 

" L 2 + + 2 ^^-^^ 

where, the effective trapping potential is defined as: 

K//(r)=yt™p(r)-^. (2.2) 

This form of the energy suggests the definition of a new norm to equip the domain 
of the functional such that the functional is coercive with respect to this norm. This 
implies that if the size of the argument is large, then naturally the value of the 
functional will be large as well, making it suitable for rotating cases. 

2.1. Inner products and norms. We define three inner products on C^(I?, C) 
and study the completion of this space with respect to the norm arising from each of 
these inner products. Consider the inner products: 

(u,w)l2 = / (2.3) 
Jv 



{u,v)h = f {u,v) + (Wu^Wv), (2.4) 
Jv 

and 

{u,v)ha^ / {u,v) + {Vau,Vav), (2.5) 
Jv 

where Va = V + ifi^*, f2 is a fixed positive number. Here (•, •) denotes the complex 
inner product. Each of these inner products leads to a norm which we will denote 
by II • IIl^JI • II H, and || • \\ha- For X ~ L'^,H,Ha, consider the completion of {u S 
C^(I?, C) : ||u||x < oo} with respect to each of the respective norms. In the first case, 
one obtains the Hilbert space = L^(I?, C), in the second case ~ H^-^{'D,IC), 
and in the third case we call the resulting Hilbert space Ha = Ha{'D,C) (see [2] 
for details on Sobolev spaces). Furthermore, the following calculation shows how the 
three norms are related. We first note that: 

{Vau,Vav) = {Vu,Vv) +Q^r^{u,v) +in{{A^u,Vv) - {Vu,A*v)). (2.6) 

If r-p denotes the radius of D, one has 

{Vau, Vau) = |Vu + inA^l^ < 2(|Vup + rl,n'^\u\'^), (2.7) 

and, consequently, 

II^^IIh^ < / (1 + 2rln')\u\^ + 2|Vu|2 < c\\u\\l (2.8) 
Jv 

where c = max{l + 2r^il^,2). Hence one has that the norm dominates the Ha 
norm. 
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In the same time, using the identity 



/ (A*u,Vw) = - / (Vw,A*t>), (2.9) 
Jv Jv 



we 



infer from (|2.6p that 



(Vau,Vaw)= / (Vm,Vu)+ / - 2il7 / (Vw,A*w). (2.10) 

Jv Jv Jv 

and, consequently, 

/ \\/u + inA^uf = / \Vuf + n^r^\u\^ -2n{iVu,A^u). (2.11) 
Also, for e > 0, one has the inequality 

ab^%e<U(^)\{be)A. (2.12) 



Now, using the above inequality and Cauchy- Schwartz one has that 

2\{i^u,A'u)\ < 2\Vu\\A'u\ < (e|Vu|)2 + ^i^. (2.13) 



Thus 



2Q{i^u,A'u) > -n{{e\Vu\)^ + ^^^) > -n{{e\Vu\)^ + ^Im^). (2.14) 



From this one has that 



/ \Vu + inA*u\^> [ |Vu|2+OV2|M|2-ri((e|VM|)2 + ^|u|2) = (2.15) 
Jv Jv £ 

/ (1 - f7e2)|VMp + (riV - n^)\u\'' > (2.16) 
Jv ^ 

[ (l-r!e2)|Vu|2-r!^|u|2 (2.17) 
Jv £ 

Now we choose e so that < 1 — Qe'^ < 1 and let fc = 1 + Since fc > 1, we can 

write 

k[ |uP + |V.uP> / k\u\^ + W Au\^ > 



f \uf + {l-ne^)\Vu\^ > (l-ne^) f \u\^ + \Vuf. (2.18) 
Jv Jv 

and infer that the Ha norm dominates the norm. Hence the two norms are 
equivalent. Furthermore, we have the following relationship between the three Hilbert 
spaces, 

H^'^{V, C) = Ha{V, C) c L^{V, C). (2.19) 

As sets and Ha are equal. However, by using the equivalent norm induced on 
Ha, we will see that the numerical performance of the descent method is improved 
for the minimization of the GP energy with rotation. 
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2.2. Gradients. The next step in writing a descent method to directly minimize 
the energy, as given in equation (|1.2p or (|2.ip . is to obtain a gradient corresponding 
to each inner product. Taking the Frechet derivative of (|1.2p . one gets that 

E'{u)h= I 3? ((Vu, V/i) + (2Ft™p w + 2g|upw- 2il7AVM,/i)) (2.20) 

JV 

or equivalently 

E'{u)h^ I {{V AU,"^ Ah) + {2Veff u + 2g\u\^u,h)). (2.21) 
Jv 

Since E'(u) is a continuous hnear functional from to R then for each u e iJ^, 
there exists a unique member of which we denote by V hE{u) so that 

E'{u)h = ^{VHEiu),h) H, (2.22) 

for all h £ ff^. We say that V hE : is a gradient for E taken with respect 

to the inner product. Likewise E'{u) is a continuous linear functional from Ha to 
R, thus it has a representation like the one given above. We denote this gradient by 
V HaE '■ Ha — >■ Ha (see [57] for a background on gradients obtained in this manner). 
Furthermore, we note from (j2.20p that for all h e 0^(1), C) one has that 

E'{u)h^^{WxE,h)x ^ I n{-\/^u + 2VtrapU + 2g\u\'^u-2inA*Vu,h). (2.23) 

Jv 

When X = L^, we directly obtain the expression of V^aiJ, the (or ordinary) 
gradient of E, already recalled in (|1.4p . From a practical point of view, it is interesting 
to note that H^ and Ha gradients will be computed using different forms of (|2.23p : 
the corresponding strong formulation for the finite difference implementation (see also 
[T7] ) and the weak formulation for the finite element implementation (see also [H]). 

3. Constrained energy minimization. 

3.1. Projection method for the mass constraint. Before discussing the 
gradient descent method, we give a brief description of the projection used to deal 
with the mass constraint. In approximating stationary states, one could in principle 
use a normalized gradient flow in conjunction with a traditional Lagrange multiplier 
for the constraint [T71III]. We adopt here a different approach and develop a projection 
method that will, in the continuous case, enforce the constraint for all time. 

The method for enforcing the constraint is presented in |27| for any general con- 
straint and hence does not provide the needed expression for our case. For the unitary 
norm constraint, several projected gradient methods are developed in [28l[7], based 
on the idea to directly compute the gradient in the tangent space to the unit sphere. 
In this work, in order to facilitate the numerical implementation, we first compute the 
gradient and then project it into the tangent space. For this purpose, it is very helpful 
to derive an explicit expression of the projected gradient that allows to preserve the 
unitary norm of the solution through the minimization procedure. It should be noted 
that explicit expressions of the projected gradient are given in [7] for the IR" gradi- 
ent flow of the linear eigenvalue problem on the unit sphere and for the gradient 
flow of the Hartree-Fock nonlinear eigenvalue problem. We derive below an explicit 
expression of the projected gradient that allows to minimize trajectories when other 
Hilbert spaces than are considered. 
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Let X = L? , H^, or Hj^. As previously stated, for each u £ X, one can find a 
member of X, denoted by \/xE{u), so that E'{u)h = {h, WxE{u))x- We called such 
an element oi X a. gradient of £' at w. Consider (5 : X ^ R, 

P{u) = / \u\\ (3.1) 
Jv 

Since we want to minimize the energy E{u) subject to the constraint P{u) = 1, we 
obtain the tangent space for our problem: 

Tu,x = null{P'{u)) = {weX : {u, w) l2 = 0}. (3.2) 

Note that T^^x is a closed linear subspacc of X , and for each u £ X , there exists a 
unique orthogonal projection from X onto T^^x- We denote this projection by Pu,x- 
Note also that Pu,x is a linear transformation with domain X and range T^^x- Thus, 
this transformation depends on the Hilbert space and u £ X. 

Let uq € X so that (3(uq) = 1 and write the descent method with the projected 
gradient: 

z(0) = uo and z'{t) = -P.^t^^x^ xE{z{t)). (3.3) 

We can easily see that /3(z) is constant since 

imnt) = f3'{z{t))z'{t) = ~p'{z{t)){P,^^)^xyxE{z{t))) = 0, (3.4) 

for all i, as Pz{t).x is the projection of X onto the null space of P'{z{t)). Thus f3{z) is 
constant and il u ~ limt^.oo z{t), then /3(m) ~ P{uq), and the norm of the initial state 
is preserved. In conclusion, by projecting the Sobolev gradient of E at z{t) into the 
null space of P'{z{t)) for each t, we get that z{t) satisfies the mass constraint for all t 
(see [27] for a more detailed development on this topic). 

For numerical implementation purposes, we give below an heuristic derivation of 
the explicit expression of the projection (see [23] for a more rigorous demonstration) . 
If, for the sake of simplicity, G = xE{u) denotes the Sobolev gradient gradient of 
E at u, the projected gradient is determined from the following two conditions: 

Pu,xQ € Tu,x, (3.5) 



{Pu,xG, h)x - E'{u)h, yh e Tu,x- (3.6) 

In order to satisfy p.6p . we choose the projected gradient of the form Pu,xG = 
Q — Bvx, with _B e R a constant and vx ^ X such as 

{vx,h)x = {u,h)L2,yheX. (3.7) 

The constant B is then obtained by imposing p.5p . The final expression that will be 
used for numerical implementation is: 

Pu.xy = y - ^ — vx, (3.8) 

k{u,vx)l^ 



with vx computed from (|3.7p . Note that if X — L^, vx = u and we recover the explicit 
expression of the projected gradient given in [6]. It is also important to note that, in 
regard to numerical consideration as well as obtaining global existence, uniqueness, 
and asymptotic convergence, we need that the map u — > Pu.x^ xE(u) be as a 
map from X to X . Using the above expression for the projection, we present in the 
next section some convergence results. 
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3.2. Convergence results in an infinite dimensional Hilbert space. In 

this subsection wc define the evolution equation we use in the Hilbert space setting 
and give our global existence and convergence result for the constrained minimization 
problem. Here we extend the results obtained in [23] for the case of the Gross- 
Pitaevskii energy without rotation. In this work, as well as in |23j . we move away 
from the general theory of Sobolev gradients as presented in [27], since the criteria 
for asymptotic convergence of the evolution equation for constrained minimization 
problems is not available in [?7] . 

The idea below the following analysis is to show that the GP energy functional 
with rotation has the same properties as the GP energy without rotation if the norm 
II • \\ha is used. We thus can adapt the results obtained in [23] to our case. We start 
by noting that, due to the mass conservation, one can add a multiple of /j, |up to the 
Gross-Pitaevskii energy and the resulting functional will have the same minimizers as 
the original functional. The idea is to obtain a functional that is uniformly and strictly 
convex. We remind the reader that for X a Hilbert space, we say that E' : X — >■ R is 
uniformly and strictly convex if there exists e > so that E"{u){h, h) > e|/i|^ for all 



Indeed, let us consider the form (|2.ip of the energy functional, and suppose that 
there exists 1 > (5 > so that K// > S. We observe that 



and infer that E : Ha — > R is uniformly and strictly convex with the assumption 
that Veff is bounded away from zero. Due to the equivalence of norms, E : H 
is also uniformly and strictly convex. Note that if Veff is not bounded away from 
zero, then one can obtain this property by adding a multiple of the constraint to the 
energy. This does not change the minimization problem as indicated by the following 
theorem. 

Theorem 3.1. Let E be a function on a subspace X contained in L^(T>). Let 



Then for (3{u) = \u\^ and h e null{(i' {u)) , E'{u)h = iff E'^{u)h = 0. 

Some other properties of the functional are required to obtain the asymptotic 
convergence of the evolution equation. In particular, we need the functional to be 
continuously twice Frechet differentiable and bounded from below. The latter two 
properties are standard and we therefore omit them. With these properties checked, 
we can give our global existence and convergence in the asymptotic limit. Using the 
space Ha, the proof of the following two theorems are identical to the ones given in 
[23j . Thus we omit the proofs and refer the reader to this work. 

Theorem 3.2. Suppose X is a Hilbert space and that _E : X — > R is continuously 
twice Frechet differentiable. Suppose also that (3 : X 9. is a given function such 
that, if Pu.x denotes the orthogonal projection of X onto the nullspace of/3'{u), then 
the map u — > Pu^x is . Then zit) given by p.3p is uniquely defined for all i > 0. 

Theorem 3.3. Suppose the hypothesis of Theorem \S.2\ and that z{t) is given by 
equation (|3.3p . with \/xE{uq) 7^ 0. If E : X ^ R is uniformly and strictly convex, 



h e X. 





(3.11) 
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then 

lim z(t) = u (3.12) 

t^oo 

exists. Furthermore, there exist two constants m and c so that \\u — z(t)\\x < me~'^*, 
and E'{u)h = for all h G null{P' {u)) . 

From the above two theorems, if we have in mind that the functional E defined 
in (|1.2p is continuously twice Frechet differentiable and uniformly and strictly convex 
when the domain is considered to be Ha or H , we obtain the result that E has a min- 
imizcr in Ha and in H^ that satisfies the constraint j3. Furthermore, this minimizer 
is obtained as the limit of the trajectory we defined in p.3p . This convergence result 
is not only important on its own, but, as we shall see, plays an important role in the 
rate of convergence of our numerical simulations. 

4. Numerical implementation. In this section we explain in detail the setup 
for our simulations using the descent method with both finite differences and finite 
elements discretization in two space dimensions. Both implementations follow the 
general lines of the algorithm described below. 

If y = Lg, Hq^ Hac denotes the finite dimensional Hilbcrt spaces resulting after 
the discretization of the domain V, the descent method p.3p takes the following 
discrete form: 

starting from uq £ Y , define a trajectory z„, n > 1 as (forward Euler scheme): 

ZO = "0, Zn+l =■ Z-a - StnVz^^^EaiZn), (4.1) 

where Vu,YEGiu) denotes the gradient obtained with respect to each inner product 
and projected following p.8|) . The time-step value Stn could be optimized when 
computed as the local minimum of the real valued function 

r^EG{z^-rW,„,YEG(zn)). (4.2) 

As convergence criterion, the algorithm stops when the relative change in energy Eg 
is below of an imposed limit. We note that in the continuous steepest descent 
algorithm, the constraint was satisfied for all time t and hence for the converged 
solution. In the discrete case, due to the first order discretization in time, it is easy 
to see from (j4.ip that the norm is conserved at time level (n + 1) up to an error 
of order ((5t„)^||Vz„^yi?G(z„)|ji2. After each (or several) iteration(s), one could also 
normalize the solution, as in j28| where a Sobolev descent method with step-size 1 
is used. This results in an improvement in the accuracy to which the constraint is 
preserved. The main observation that we made was that even though we used a first 
order discretization in time, our projection method allowed to take larger time steps 
when compared to the method using the normalization alone. 
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4.1. Finite differences. Wc discretize T) into an by equally spaced [S^ = 
5y — 5) grid and let Vq be the set of all K — N"^ grid points. Let X be the collection of 
all complex valued functions on Dq. For / G X, (Dif)(x,y) is the approximation to 
the partial derivative in the first independent variable at (a;, y) and {D2f){x, y) is the 
approximation to the partial derivative in the second independent variable at {x,y). 
We have used a fourth-order centered finite-difference scheme to approximate the first 
partial derivatives. When compared to the classical second order scheme, this high- 
order approximation proved very helpful in computing complex configurations (vortex 
lattices within the condensate) with reasonably fine grids. Furthermore Df = (^^z)' 
For {x, y) a grid point, we also define Di^a and D2,a by 

(i?i,^/)(x, y) - (i?i/)(x, y) + inyf{x, y) (4.3) 

and 

{D2M)[x^ y) = {D2f){x, y) - tflxfix, y). (4.4) 

We denote by Da = (^j'^); the discretized form of the operator V a- The three inner 
products that equip X are defined as: for f,g d X, 

{f,9)L^^{f,9), (4.5) 



{f,9)H = + {Dif,D,g) + {D2f,D2g), (4.6) 

and 

(/, .9)^.4 = if, a) + {D,,Af,D,,Ag) + {D2,Af,D2,A9), (4.7) 

where (•, •) denotes the complex inner product. Note that (•, ■)l2 is analogous to 
the L^(I?, C) inner product, (., ■)h is analogous to the C) inner product, and 

(., ■)ha is analogous to the i/^(2?, C) inner product. 

Since Di, D2, Di a, D2^a can be viewed as a linear transformation acting on 
C^^, we think of each of these transformations &s &. K x K matrix. Let D\[ denote 
the conjugate transpose of the corresponding matrix. Wc note that we can write the 
H and Ha inner products as 

(/, g)H = ((/ + DID, + D*Mf.9)L^- (4.8) 

and 

= {{I + Dl^D,.A + DI^D2,a)L9)l-- (4.9) 

The collection X makes a finite dimensional Hilbcrt space with each of the above 
inner products. We denote the resulting Hilbcrt spaces by L^, Hq, and Haq- Now, 
we discretize the energy functional as given in equations (|1.2p and (|2.ip . Here the 
subscript G denotes that we are in the finite difference setting. 

Ecif) = 'J'E ^d^i/l' + + Vtrapo\f? + f I/I' -^rotaf. (4.10) 

■Dg 

where for x E Vg, 

rota fix) = (*/(a;)M(rr) [^j^'J^^^^^) , (4-11) 
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and 

A{x,y)^{y -x). (4.12) 

Equivalently, 

Eaif) ^5^T.\ {\Di,aI? + \D2A?) + ^e//, |/P + fl/l^ (4.13) 

■Dg 

If we take a derivative of Eq, we get that 

■Dg 

+ 2(/i, + - f}a*(D/)). (4.14) 

Observe that for each f G X, E'Q{f) is a continuous hnear transformation on X using 
any of the three norms we specified. Thus it has a representation with respect to each 
of the inner products we defined above. Using this representation, we wiU obtain a 
gradient. Since the inner product is proportional to the Eucfidean inner product, 
the ordinary or Eucfidean gradient (i.e. the hst of partial derivatives of Eq taken with 
respect to each of the K independent variables) is easily derived if the real valued 
transformation E'Q{f) is rewritten as: 

E'a{f)h^^{h,^L-E{J))L^. (4.15) 

We get that 

V L-E{f) ^ 8\DlDif + D;D2f +2iVtrapGf +9\f\^f - niA\Df))).{4.16) 

We now derive the other two gradients, V HE{f) and S/HAE{f), with respect to 
the H and Ha inner products. From (|4.8p and (|4.9p we obtain that 

E'a{f)h = ^{h, VHEG{f))H = "^{h, {I + D*D)VHEG{f))L-, (4.17) 

and 

E'a{J)h^^{h,V HAEG{f))HA ^"^{^{1 + D\Da)V HaEgU))l- (4.18) 
By comparing these equations to (|4.15p . we finally get that 

"^hEgU) = {I + D*D)-'VL2E{f), (4.19) 

and 

"^H^EGif) = (/ + D^DAr'^^L-Eif). (4.20) 

The discrete descent method (|4.ip using the above finite difference dscretization 
was implemented in Matlab. The Sobolev gradients are computed from (|4.19p and 
()4.20p by solving linear systems at each time step using a preconditioned conjugate 
gradient method. Since this part is time consuming on fine grids, we used a linesearch 
algorithm to locally compute the time step from (j4.2p . This resulted in a significant 
reduction of the number of iterations needed to achieve convergence. 
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4.2. Finite elements. The finite-clcmcnts implementation uses the free soft- 
ware PreeFem-l— I- [5D], whieh proposes a large variety of triangular finite elements 
(linear and quadratic Lagrangian elements, discontinuous PI, Raviart-Thomas ele- 
ments, etc.) to solve partial differential equations (PDE) in two dimensions (2D). 
FreeFemH — h is an integrated product with its own high level programming language 
with a syntax close to mathematical formulations. 

It is therefore very easy to implement the variational formulations associated to 
the calculation of the three gradients, since the definitions of scalar products (|2.3p - 
(|2.5p use an integral form. Following the developments in section (|2.2p . and using 
as definition of the complex inner product {u,v) = uv* , the ordinary gradient is 
derived from (|2.20[) and computed as the solution Q = '^l^E of the problem with 
homogeneous Dirichlet boundary conditions: 

/ gh = miS, (4.21) 
Jv 

RHS= / WuWh + 2[VtrapU+ {g\u\^)u-inA^Wu]h, (4.22) 
Jv 

where h stands now for the real valued basis function of the finite element space. 
Following (|2.23p . the gradient is directly computed by solving the equation: 

/ VgWh + gh^ RHS, where G = VhE. (4.23) 
Jv 

It is interesting to note that (|4.23p is directly derived from the weak formulation of 
(|2.23p , with the obvious advantage to obtain a simpler right-hand side (|4.22p . which is 
derived by integrating by parts the weak form of the gradient. Therefore, in order 
to solve (|4.23p . it is not necessary to explicitly compute the gradient (by solving 
(|4.2ip ). as required for the finite-difference implementation. 

Observing from (|2.10p that the Ha scalar product could be expanded to obtain 
the equivalent definition: 

<u,w>H^= / {[l + n^iy^ + x^)\u,v) + {\/u,\/v) ~2in{A^\/u,v), (4.24) 
Jv 

the Ha gradient is directly computed as the solution Q = V HaE of the problem: 

/ [i + n^iy^ + x^)\gh + \/gvh-2in{A^\/g)h = Rm. (4.25) 

Jv 

It is interesting to emphasize the fact that previous equations are solved in complex 
variables. The approach based on the separation of the real and imaginary part of the 
gradient used in [29] is not possible when computing the Ha gradient. The FreeFem 
scripts are written in an optimized form using the pre-computation and factorization 
of the complex matrices associated to linear systems given by (|4.23p and (|4.25p . It is 
interesting to note that the same matrices are involved in the computation of vx from 
(|3.7p : the projected gradient (|3.8p could be therefore optimized in the same way. The 
implementation uses PI (piecewise linear) finite-elements, with a P4 representation 
of the nonlinear terms appearing in (j4.22p . A fifth order quadrature formula was 
used to compute two-dimensional integrals. The FreeFem scripts allow to switch to 
P2 (piecewise quadratic) finite elements by a simple change of the definition of the 
generic finite-elements space. Adaptive mesh refinement was used for simulations of 
rotating BEG with dense lattice of vortices. 
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5. Numerical experiments. Wc first use a test case with analytical manufac- 
tured solution to ascertain the convergence of the steepest descent method for each of 
the three gradients. Then, we use the numerical set-up to compute simple metastable 
states of rotating Bose-Einstein condensates with single or multiple vortices. The 
performances of the three methods are comparatively evaluated. Finally, the new Ha 
gradient method is used to compute complex configurations relevant for real rotating 
condensates (Abrikosov vortex lattice and giant vortex). 

5.1. Test case with manufactured solutions. This test case is used as bench- 
mark for the evaluation of the descent method for each of the three gradients (L^, H, Ha)- 
Wc consider a non-linear problem close to the Gross-Pitacvskii equation: 

- iv^u + Ctrap u + g\u\^u - in{A'\/)u = /, (5.1) 
corresponding to the minimization of the energy functional: 

E{u, f)^ ^|Vu|2 + Ctrap \uf + 5 ^ - + ful - n^{iu*A*Vu). (5.2) 

For this energy functional, the gradient is expressed as in (j2.23p . with Vtrap = 
Ctrap = const, and a supplementary term — 2(/, h) to be added. It should be noted 
that this is a test case of minimization without constraint. 

In order to test the implemented methods, we manufacture solutions of (|5.1|) : 
we consider a given expression for u and calculate the corresponding right-hand side 
f{x, y). A simple way to construct such manufactured solutions is to consider solutions 
with azimuthal symmetry: 

Uf{x,y)~U{r)c-K.p(im9), (5-3) 



where (r, 6*) are cylindrical coordinates (r — \J + y"^). Since the Laplacian in cylin- 
drical coordinates reads 

V2-lAfr— V— — (5 4) 

r dr \ dr ) 89"^ ' 

and the new term corresponding to the rotation becomes 

du du du . ^. 



we obtain that 



/ = F(r) exp(im6l), (5.6) 



with 



= -^-T- f^?^) +1"^^ + ^trap U + gU'~ mW. (5.7) 
2 r or \ or I 2 



Wc choose the domain I? to be a circle of radius R and 



U ^ r^{R^r), 



(5.8) 
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which satisfies the homogeneous boundary condition u = for r — R. For this choice, 
we obtain useful analytical formulas for 

F[r) = -i(4i? - 9r) + ]^m^{R - r) + Ctrap U + CnU^ - mW, (5.9) 
and energy 

/ 7->6 d6 d8 "^ri^ \ d8 

E(u, /) = 27r -— -m^—~ Ctrav— - Cn-— + mVLn—. (5.10) 
^ V 20 120 ^168 20020y 84 ^ ' 

The contour patterns for such solutions are displayed in Fig. l5.1l for m — 1 and m = 3. 



m=1,Real(u) ^ m=1, Imag(u) 




Fig. 5.1. Contour patterns of the manufactured solution corresponding to equations 1)5. 3p and 
1 15.811 . Azimuthal wave numbers m = I and m = 3. 

The numerical application for manufactured solutions consider the following pa- 
rameters: 

Ctrap = 20, g = 100, R = l, m = 3, = 10. 

For this case, the theoretical values for energy and angular momentum of the exact 
solution are: E = —0.505553 and = 0.1122, respectively. The computation is 
considered as converged if the relative variation of the energy is less than e = 10^^. 

Tables 15.11 and 15.21 assess the convergence of the descent method by computing 
different norms of the difference between the exact and computed solutions. Perfor- 
mance of each gradient method are quantified by extracting the overall computing 
(CPU) time and the number n of time steps necessary to achieve convergence. 

All test cases considered mq = as the initial guess for the descent method. 
Different initial conditions (e.g. uq computed as the solution of the corresponding 
linear problem) were tested with similar convergence results. 
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N 


n 


CPU 


E(u) 


1 M ~ ^i/lloo 


I kt — M f 1 r 2 

II J 1 1 -tj 


I M — Uf llff 

II J 1 1 ^-^ 




2^ 


1137 


35.41 


-.4828 


.0010 


2.25e-4 


.0270 


H 




100 


16.23 


-.4828 


2.81-4 


1.66e-5 


.0021 


Ha 


2« 


38 


6.39 


-.4828 


1.41C-4 


3.85C-6 


4.73e-4 


L 


2' 


1960 


308.13 


-.4941 


4.11C-4 


2.84e-4 


.0185 


H 


2' 


40 


36.34 


-.4941 


6.32e-5 


4.81e-4 


5.28e-4 


Ha 


2> 


18 


17.22 


-.4941 


2.05e-5 


6.84e-7 


6.26e-5 


L 


2« 


> 3000 


> 2.04e3 


-.4985 








H 


2« 


30 


154.47 


-.4997 


4.59-5 


9.77e-6 


.0013 


Ha 


2« 


14 


73.21 


-.4997 


1.56G-6 


1.36e-6 


1.504e-4 



Table 5.1 

Test case with manufactured solutions. Algorithm efficiency and convergence test for the finite 
difference implementation (variable time step computation). 
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CPU 


E(u) 


\\u - U/ oo 


\\u-Uf\\L2 


\\u-Uf\\H 


5t 


L 


100/1776 


1176 


85 


-.4934 


1.988C-3 


l.OOlc-3 


1.705e-2 


8e-4 


H 


100/1776 


47 


3.4 


-.4934 


1.883e-3 


9.220e-4 


1.668e-3 


1 


Ha 


100/1776 


14 


1 


-.4934 


1.880e-3 


9.140e-4 


1.665e-2 


3 


L 


200/7064 


4292 


1252 


-.5025 


7.492e-4 


4.200e-4 


7.401e-3 


2c-4 


H 


200/7064 


47 


13.8 


-.5025 


5.530e-4 


2.232e-4 


6.548e-3 


1 


Ha 


200/7064 


14 


4.1 


-.5025 


5.390e-4 


2.119e-4 


6.474e-3 


3 


L 


400/27604 


> 8000 


> 9193 


-.5027 








5c-5 


H 


400/27604 


47 


54.2 


-.5047 


1.687e-4 


6.8535e-5 


3.954C-3 


1 


Ha 


400/27604 


14 


16.2 


-.5047 


1.549e-4 


5.730C-5 


3.791e-3 


3 



Table 5.2 

Test case with manufactured solutions. Algorithm efficiency and convergence test for the finite 
element implementation (fixed time step computation). The triangular mesh is generated with M 
points on the border of the domain. 



The first obvious observation is that the descent method using the ordinary 
gradient has very slow convergence rate because of very small time steps imposed by 
the stability limit of the method. This was expected since this method is the equivalent 
to the explicit Euler integration scheme for the imaginary-time propagation equation. 
A similar result was reported in |29j for simpler Schrodinger type energy functionals. 
Larger time steps are allowed in the and Ha methods, since the Sobolev gradients 
represent a preconditioning of the ordinary gradient |171 [^t i6j . 

For the descent methods using a constant time step 5t (finite element implemen- 
tation), we compare the computations performed using the maximum value {5t)max 
allowed by the stability of each method. These values, displayed in Tab. 15.21 were ob- 
tained by successive tests: the value of 6t was increased by 20% for each new run, until 
the computation became unstable. It should be noted that we were not interested in a 
refined numerical evaluation of the stability limit of each method, since computations 
using a more precise estimation of {St)max did not result in a significant variation of 
the CPU time. The same approach to compare methods using their maximum time 
step allowed by stability reasons will be applied to all subsequent computations in 
this section. 

Tables 15.11 and 15.21 also allows to relate the computing cost to the complexity of 
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each method. As already stated, the descent method using the gradient can be 
regarded as an exphcit backward Euler scheme. It therefore has httle complexity and 
the computing cost per iteration step (i.e. the ratio CPU/n) is very low. Sobolev 
gradients are computing by solving linear systems, which adds extra computational 
cost. For the finite-difference implementation, equations (|4.19[) or (|4.20[) are solved by 
a preconditioned conjugate-gradient method; since this part of the algorithm is time 
consuming, the CPU time per iteration step (CPU/n) is multiplied up to a factor 
of 8, when compared to the gradient method. The situation is different in the 
finite-element implementation. Since the weak formulation of the equation (|2.23p is 
used, the computation of all gradients needs to solve a linear system. In order to have 
an optimized numerical implementation that can switch between the three descent 
methods, the matrix of this system is stored and factorized before the time loop. As a 
consequence, even though the matrix of the system in (|4.2ip is simpler (mass matrix) 
than in (|4.23p or (|4.25p . the ratio (CPU/n) is identical for the computation off all 
gradients. 

In all numerical tests, the convergence of the gradient method needs a large 
number of time steps, and, consequently, much larger CPU times than the Sobolev 
gradient methods. Since the performances of the gradient method are very poor, 
it will not be used in the following numerical experiments. We shall now focus on the 
comparison between the H and Ha method. For this test case considering a large value 
of f2, the Ha gradient allows for larger time steps and therefore the computational 
time is considerably reduced, by approximately a factor of 3. This suggests that 
the preconditioning of the gradient introduced by the new Ha inner product is very 
effective for computing eases with high rotation frequencies VL (it goes without saying 
that the H and Ha methods are equivalent for — )■ 0). 

5.2. Simulations of rotating Bose-Einstein condensates. In computing 
stationary states of rotating Bose-Einstein condensates, the initial state uq in the 
descent method (|4.1[) plays a crucial role. The algorithm usually starts from a wave 
function distribution derived from the Thomas-Fermi approximation. In the strong 
interaction regime (large values of g), it is reasonable to neglect the contribution of 
the kinetic energy and work with the simplified energy functional: 



The minimizer of this energy corresponds to the Thomas-Fermi atomic density: 



where /i is the chemical potential. Since /i is a Lagrange multiplier, imposing the 
mass constraint in (j5.12p yields a relation for /i. After computing the value of /i, the 
Thomas- Fermi radius of the condensate can be determined from (|5.12p {ptf{Rtf) = 
0). When a rotation 51 is applied, the Thomas- Fermi approximation (j5.12p stands 
with Ve// replacing Vtrap- The resulting radius R^p is used to estimate the size of 
the domain T) in simulations {ro > Rtf) ■ 

Wc also mention that the converged final state is characterized by its energy E{u) 
and angular momentum Lz{u) which gives a measure of the rotation: 




(5.11) 




(5.12) 




(5.13) 
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5.2.1. Off-center vortex case: harmonic trapping potential and small 

n. The second numerical experiment considers the classical harmonic trapping po- 
tential and an initial state computed from the Thomas-Fermi approximation plus a 
singly quantized vortex of center located at (x^, Uv)- We use an ansatz for the vortex 
described in [3]. The parameters of the simulation are the following: 



g = 500, Vt 



trap 



72, 17 = 0.4, 



0.5, y„ = 0. 



(5.14) 



The Thomas-Fermi radius is for this case R^p 
is circular of radius R ~ \.2bR^p 



5.246 and the computational domain 
6.56. The final converged state contains a single 
vortex centered at the origin (see Fig. 15. 




o 



Fig. 5.2. Off-center vortex case. Initial state with an off-center vortex and final converged state 
with a centered vortex. Contours of atomic density 
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E(u) 


LM 
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2« 
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169.51 


8.3587 


.9998 


Ha 


2« 


1197 


166.34 


8.3587 


.9998 


H 


i' 


1184 


866.88 


8.3605 


.9999 


Ha 


i' 


1127 


890.06 


8.3605 


.9999 


H 


2« 


1274 


4.9548e3 


8.3606 


.9999 


Ha 


2« 


1244 


4.7882e3 


8.3606 


.9999 



Table .5.3 

Off-center vortex case. Algorithm efficiency and characterization (E{u), L^{u)) of the converged 
state state for the finite difference implementation (variable time step computation). 



The comparative results are presented in Tabs. 15.31 and 15.41 It is important to 
note that the convergence test must be set to e = 10^^ in order to obtain a final 
state with a vortex centered at the origin and = 1 (theoretical value reached for 
the finest meshes). A relaxed convergence criterion will result in a vortex that is not 
exactly centered since the convergence rate is very slow at the end of the simulation. 
As expected, the Hi and Ha perform similarly because of the low value of il. 



18 



I. DANAILA AND P. KAZEMI 





M/triangles 
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CPU 


E(u) 




H 


100/1762 


701 


56.97 


8.3819 


.994598 


Ha 


100/1762 


703 


57.47 


8.3795 


.994575 


H 


200/7064 


1667 


537.85 


8.3720 


1.00042 


Ha 


200/7064 


1717 


556.41 


8.3694 


1.00022 


H 


400/27604 


1788 


2.335c3 


8.3699 


1.00052 


Ha 


400/27604 


1831 


2.407e3 


8.3673 


1.00032 



Table 5.4 

Off-center vortex case. Algorithm efficiency and characterization (E(u), Lz (u) ) of the converged 
state state for the finite element implementation. The time step is set to 0. 1 for all computations. 



5.2.2. Vortex array case: harmonic-plus-quartic trapping potential and 
large fi. The harmonic trapping potential physically sets an upper bound for the 
rotation frequency, since for = 1 the centrifugal force balances the trapping force 
and the confinement of the condensate vanishes. To overcome this limitation, different 
forms of the trapping potential are currently experimentally and theoretically studied. 
We use in the third numerical experiment a combined harmonic-plus-quartic potential 
(see also [21] |4j [131 [IS] ) with the following parameters 



5 = 500, + Q 



(5.15) 



The Thomas-Fermi radius is for this case 



3.40. The computational domain is 
circular of radius Rmax = 1.25i?^^. The initial state contains a central vortex plus 
an array of 6 vortices equally distributed on the circle of radius 0.25i?,„a^. All the 
vortices have a winding number m = 1, except the first vortex that has m = 2 (Fig. 
5.3[) . Since vortices with winding number m > \ are not physically stable, the m = 2 
vortex will split into two singly quantized vortices. The final state contains therefore 
a central vortex an array of 7 vortices (Fig. 15. 3p . The convergence test is relaxed to 
e = 10-^ 



-41 

4 -4 



Fig. 5.3. Vortex array case. Initial state with 6 vortices and final converged state with an array 
of 7 vortices. Contours of atomic density \u\'^ . 



Tables 15.51 and 15.61 show that the converged state is the same for both finite 
difference and finite element implementations. The Ha method has better stability 
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gradient 
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2.245e3 
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6.4691 



Table 5.-5 

Vortex array case. Algorithm efficiency and characterization (E{u), L:;{u)) of the converged 
state state for the finite difference implementation (variable time step computation). 



gradient 


Af 1 triangles 


n 


CPU 


E(u) 




H 


100/1762 


507 


42.28 


12.0553 


6.1297 


Ha 


100/1762 


330 


27.70 


12.1413 


6.1654 


H 


200/7064 


418 


138.53 


11.5341 


6.3920 


Ha 


200/7064 


270 


90.11 


11.6171 


6.4135 


H 


400/27604 


420 


550.10 


11.4017 


6.4641 


Ha 


400/27604 


262 


346.87 


11.4846 


6.4840 



Table 5.6 

Vortex array case. Algorithm efficiency and characterization (E{u), Lz(u)) of the converged 
state for the finite element implementation. The maximum allowed time step is 0.1 for the H 
gradient and 0.2 for the Ha gradient. 



properties and allows a CPU time gain up to 36%. This gain was expeeted sinec 17 is 
large for this case. 

5.2.3. Giant vortex and Abrikosov vortex lattice. Finally, to show that 
the new method has the capability to handle more complicated cases, we produce the 
giant vortex using the Ha gradient in conjunction with the new projection method 
proposed to enforce the mass constraint. We use the parameters of the previous 
numerical experiment (harmonic-plus-quartic trapping potential) and progressively 
increase 51 from 2 to 4. Each computation starts from an initial field representing 
the converged state previously obtained for a lower value of VL. The transition from a 
vortex lattice to the giant vortex is observed (Fig. 15. 4p . The giant vortex is a hole in 
the condensate (the atomic density goes to zero inside) with multiple phase defects. 
This particular vortex structure, theoretically analyzed in numerous studies j211 Ul 
I13[ I15j . was captured using both the finite elements and finite difference simulations. 

A last complex computational case is illustrated in Fig. 15.51 For a harmonic 
trapping potential and high rotation frequency (fi = 0.95) an Abrikosov vortex lattice 
forms in the condensate. The difficulty in computing this case in the strong-interaction 
regime (large values of g) comes from the fact that the condensate becomes larger and 
the vortex lattice denser when the value of g is increased. In order to increase conver- 
gence, each computation starts from an initial field representing the converged state 
obtained for a lower value of g. During the iterative process, new vortices nucleate at 
the boundaries and slowly move towards their final equilibrium locations. In comput- 
ing such configurations, containing several hundreds of vortices, the adaptive mesh 
refinement capabilities of FreeFem proved very helpful in reducing the computational 
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Fig. 5.4. Giant vortex case. Converged states for = 2.5, 3, 4 showing the formation of a hole 
in the condensate (giant vortex) for high rotation rates Contours of atomic density 



time and correctly capturing vortex positions. 



g= 1000 g= 5000 g= 10000 




-15 -ID -5 5 10 15 -1 5 -1 -5 5 10 1 5 -1 5 -1 -5 5 1 15 



Fig. 5.5. Abrikosov vortex lattice case. Converged states for Q = 0.95 and increasing values 
of the interactions constant g. Finite elements computations using mesh adaptivity. Contours of 
atomic density \u\'^ . 

6. Summary. The numerical study of a rotating Bose-Einstein condensate has 
been the subject of many numerical studies, both in two dimensions (2D) and three 
dimensions (3D). Since most of the studies [SJ [HJ [HI [31 IH [131 [5J [TT] use the imaginary 
time propagation method (equivalent to the gradient flow model (|1.4p ). there are 
few studies using direct minimization by Sobolev gradient methods. Nevertheless, 
replacing the ordinary gradient in a descent method with the Sobolev gradient 
proved effective in minimizing the 3D Gross-Pitaevskii energy |181 117] or simpler 
Schrodinger type functionals [29] . 

In this work we introduced a new inner product (Ha) to equip the domain of 
the GP energy functional with rotation and derived the corresponding gradient. We 
demonstrated that numerical performance is enhanced by replacing in the descent 
method the or gradients with the gradient obtained from the Ha inner product. 
The gain in computational time proved very important when configurations with high 
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rotation rates are computed. We also introduced a new projection method to enforce 
the mass constraint. This method avoids more comphcated approaches using an 
energy functional with a penalty term, or the traditional normalization method that 
performs the descent over a path of functions with an imposed norm. 

These two new tools allowed to implement robust descent methods using finite 
difference and finite element spatial discretization. Both numerical settings proved 
very efficient in computing various complex two-dimensional configurations of rotating 
Bose-Einstein condensates. 

We finally emphasize the fact that the new gradient and projection method for 
the mass constraint have a more general interest and could be also used in conjunction 
with existing numerical schemes (such as sophisticated time stepping procedures) to 
study the energy minimization of Gross-Pitaevskii type functionals. 
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